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We  present  a  new  method  for  stability  and  modal  analysis  of  shear  flows  and  their  acous¬ 
tic  radiation.  The  Euler  equations  are  modified  and  solved  as  a  spatial  initial  value  problem 
in  which  initial  perturbations  are  specified  at  the  flow  inlet  and  propagated  downstream  by 
integration  of  the  equations.  The  modified  equations,  which  we  call  one-way  Euler  equa¬ 
tions,  differ  from  the  usual  Euler  equations  in  that  they  do  not  support  upstream  acoustic 
waves.  It  is  necessary  to  remove  these  modes  from  the  Euler  operator  because,  if  retained, 
they  cause  instability  in  the  spatial  marching  procedure.  These  modes  are  removed  using 
a  two-step  process.  First,  the  upstream  modes  are  partially  decoupled  from  the  down¬ 
stream  modes  using  a  linear  similarity  transformation.  Second,  the  error  in  the  first  step 
is  eliminated  using  a  convergent  recursive  filtering  technique.  A  previous  spatial  march¬ 
ing  method  called  the  parabolized  stability  equations  uses  numerical  damping  to  stabilize 
the  march,  but  this  has  the  unintended  consequence  of  heavily  damping  the  downstream 
acoustic  waves.  Therefore,  the  one-way  Euler  equation  could  be  used  to  obtain  improved 
accuracy  over  the  parabolized  stability  equations  as  a  low-order  model  for  noise  simulation 
of  mixing  layers  and  jets. 


I.  Introduction 

The  high  cost  of  direct  and  large-eddy  simulation  of  jet  noise  has  prompted  the  development  of  reduced- 
order  methods  based  on  instability  waves.  Experiments  have  long  suggested  that  the  breakdown  of  coherent 
structures  at  the  end  of  the  potential  core  is  the  dominant  source  of  sound  at  small  angles  to  the  jet  axis, 
where  sound  is  most  intense. 1-4  These  coherent  structures  can  be  quantitatively  identified  as  instability 
waves  in  both  forced5  and  natural6  jets.  Therefore,  accurate  modeling  of  instability  waves  has  the  potential 
to  provide  an  accurate,  low-cost  jet  noise  model. 

The  most  inexpensive  instability  wave-based  approaches  are  locally  parallel  methods  based  on  linear 
stability  theory.7,8  These  methods  have  provided  insight  into  the  role  of  the  Kelvin- Helmholtz  instability  in 
jet  noise  production  and  the  sensitivity  of  this  instability  to  various  flow  parameters.  However,  the  locally 
parallel  assumption  of  these  methods  limit  their  ability  to  accurately  capture  nonparallel  effects  and  the 
downstream  evolution  of  the  instability  waves. 

To  fully  account  for  nonparallel  effects,  global  mode  analysis  has  been  used9  to  capture  the  spatial  growth 
and  decay  of  near-field  instability  waves  in  jets  and  their  associated  acoustic  radiation.  The  trade-off  is  that 
global  modes  are  expensive  to  compute,  limiting  their  utility  as  a  reduced-order  model. 

In  between  these  two  extremes  lies  an  ad  hoc  generalization  of  linear  stability  theory  call  the  parabolized 
stability  equations  (PSE).10  PSE  can  partially  capture  nonparallel  and  nonlinear  effects,  and  has  been  shown 
to  accurately  model  many  convectively  unstable  flows.  In  particular,  our  group  has  previously  shown  that 
linear  PSE  can  produce  quantitatively  accurate  near-field11  and  reasonable  far-held12  sound  predictions  for 
unforced  turbulent  jets  at  a  fraction  of  the  cost  of  LES  or  global  mode  analysis. 

The  efficiency  of  PSE  is  achieved  by  using  a  spatial  marching  technique  in  the  streamwise  direction.  Initial 
conditions  are  specified  at  the  jet  inlet  and  propagated  downstream  by  integration  of  the  PSE  equations. 
Because  of  this  one-way  solution  strategy,  upstream  propagating  waves  cannot  be  captured  by  PSE.  This  is 
a  minor  loss,  however,  because  of  the  substantial  downstream  directivity  of  jet  noise. 

Care  must  be  taken  to  assure  a  stable  downstream  march.  This  is  the  case  because,  despite  their  name, 
the  parabolized  stability  equations  are  not,  in  fact,  parabolic.  Even  though  PSE  cannot  capture  upstream 
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sound,  the  PSE  operator  still  contains  upstream  acoustic  modes  which  constitute  an  elliptic  effect  in  the 
frequency  domain.  These  elliptic  modes  are  a  consequence  of  the  inherent  boundary  value  nature  of  the 
subsonic  Euler  equations.  Specifically,  decaying  upstream  evanescent  modes  are  wrongly  interpreted  by  the 
PSE  operator  as  growing  downstream  inodes,  causing  exponential  instability  in  the  march  if  they  are  not 
eliminated.  Eliminating  these  elliptic  modes  is  the  goal  of  parabolization.  The  traditional  PSE  method 
accomplishes  this  by  using  an  implicit  Euler  axial  discretization  along  with  a  restriction  on  the  minimum 
step  size.  This  has  the  effect  of  numerically  clamping  the  upstream  acoustic  modes  such  that  they  move 
inside  the  unit  circle  in  the  discrete  spectral  space.  If  a  step  size  below  the  restriction  is  used,  these  modes 
again  move  outside  of  the  unit  circle  and  cause  unbounded  instability.13 

An  unintended  consequence  of  the  numerical  damping  used  to  eliminate  the  upstream  acoustic  modes  is 
that  the  desired  downstream  acoustic  modes  are  also  heavily  damped.  Additionally,  the  minimum  step  size 
restriction  in  some  cases  prevents  adequate  convergence  of  the  solution.  Several  alternative  methods14, 15 
have  been  devised  to  remove  the  upstream  acoustic  mode  from  the  PSE  operator  and  reduce  the  minimum 
allowable  step  size,  but  each  of  these  methods  still  has  catastrophic  effects  on  the  downstream  acoustic 
radiation. 

In  this  paper,  we  introduce  a  new  spatial  marching  technique  that  can  be  used  to  model  instability 
waves  and  their  associated  acoustic  radiation.  Our  method  explicitly  removes  the  upstream  acoustic  modes 
from  the  linearized  Euler  equations  while  accurately  retaining  the  downstream  acoustic  modes.  The  resulting 
equations  can  be  integrated  without  the  use  of  numerical  damping,  allowing  the  downstream  acoustic  field  to 
be  accurately  modeled.  Our  method  draws  from  techniques  for  generating  nonreflecting  boundary  conditions. 
Guddati16  developed  a  one-way  marching  technique  for  the  elastic  wave  equation  based  on  a  different  class 
of  nonreflecting  boundary  conditions,  but  it  has  not,  to  our  knowledge,  been  extended  to  equations  for  fluid 
dynamics. 


II.  Method 

In  this  section,  we  outline  our  method  for  parabolizing  the  Euler  equations.  We  first  discuss  relevant 
aspects  of  the  nature  of  the  subsonic  Euler  equations,  then  present  the  two  steps  of  our  method,  and  finally 
consider  computational  cost.  We  will  focus  on  the  two-dimensional  Euler  equations  for  ease  of  presentation, 
but  our  method  can  be  directly  applied  to  three  dimensional  Cartesian  and  cylindrical  coordinates  as  well, 
and  can  be  extended  to  the  Navier-Stokes  equations. 


II. A.  The  linearized  Euler  equations 

We  begin  with  the  two-dimensional  compressible  Euler  equations: 
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The  fluid  is  a  perfect  gas  with  specific  heat  ratio  7.  We  have  chosen  as  independent  variables  spe¬ 
cific  volume,  axial  velocity,  transverse  velocity,  and  pressure.  All  variables  have  been  appropriately  non- 
dimensionalized  by  an  ambient  sound  speed  and  density  and  a  problem  dependent  length-scale.  We  decom¬ 
pose  the  flow- field  q  =  {v,  ux,  uy,p}7  into  a  steady  baseflow  and  fluctuations  about  this  baseflow: 


q{x,  y ,  t )  =  q{x,  y)  +  q'(x,  y ,  t). 


(2) 


If  non-linear  products  of  fluctuations  are  neglected  and  the  baseflow  satisfies  equation  (1),  then  the 
fluctuations  satisfy  the  linearized,  two-dimensional,  compressible  Euler  equations: 
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Next,  the  fluctuations  are  further  decomposed  into  a  summation  of  frequency  modes: 


(4) 


q'(x,y,t)  =  ^2qn{x,y)e  lUnt.  (5) 

n 

In  choosing  ojn  to  be  real-valued,  we  restrict  our  analysis  to  flows  that  are  stationary  in  time.  Inserting  this 
decomposition  into  equation  (3),  we  have  for  each  frequency  mode 

-iu>nqn  +  A(q +  B(q +  L(q)qn  =  0.  (6) 

Since  we  wish  to  interpret  this  equations  as  a  spatial  evolution  equation  in  the  streamwise  x  direction, 
we  rewrite  equation  (6)  as 


Qq 

=  -A-1  ( -iunI  +  BDy  +  L)qn  :=  Mnqn.  (7) 

The  inverse  of  the  matrix  A  exists  as  long  as  0  <  ux  <  1,  so  the  baseflow  in  the  axial  direction  must 
be  downstream  and  subsonic  everywhere.  The  operator  Dy  represents  the  transverse  derivative,  and  will  be 
approximated  by  a  differentiation  matrix  upon  discretization  of  the  equations.  The  streamwise  evolution  of 
equation  (7)  is  governed  by  the  eigenvalues  and  eigenvectors  of  the  operator  Mn.  Analytical  solutions  for 
these  quantities  do  not  exist  for  a  general  spatially  inhomogeneous  baseflow. 

To  understand  the  structure  of  Mn,  it  is  helpful  to  look  at  the  special  case  of  a  uniform  baseflow,  for 
which  analytical  expressions  can  be  derived  by  taking  a  Fourier  transform  in  y.  Then,  the  eigenvalues  and 
left  eigenvectors  are 
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(9) 


Here,  Mx  =  ux/c  is  the  Mach  number,  d  =  vjc.  is  the  acoustic  admittance,  k  =  w/c  is  the  reduced 
frequency,  c  is  the  local  sound  speed  of  the  baseflow,  q  is  the  transverse  Fourier  wavenumber,  z  =  q/k  is  a 
scaled  transverse  wavenumber,  and  the  function  y,(z)  is  given  by 
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Kz)  =  V1-  (1-  M%)z2. 


(10) 

The  first  two  eigenvalues,  As  and  Xu,  describe  the  evolution  of  entropy  and  vorticity,  respectively.  These 
modes  convect  downstream  with  group  velocity  ux.  The  eigenvalues  A-t  represent  acoustic  branches  of  the 
spectrum,  consisting  of  acoustic  waves  propagating  in  all  directions  as  well  as  evanescent  waves  which  decay 
exponentially  in  y.  It  can  be  easily  verified17  that  for  all  z,  A±  have  positive  and  negative  group  velocities, 
respectively.  These  eigenvalues  are  shown  graphically  in  figure  1.  The  parts  of  the  acoustic  branches  that 
lie  along  the  stability  boundary  are  the  propagating  acoustic  waves  and  the  vertical  parts  are  the  evanescent 
modes.  Although  the  preceding  analysis  is  strictly  valid  only  for  a  uniform  baseflow,  inhomogeneous  flows  like 
jets  and  mixing  layers  contain  modes  that  can  be  grouped  into  families  with  the  same  qualitative  properties 
as  the  four  uniform  flow  modes. 
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Figure  1:  Uniform  flow  Euler  spectrum. 


Overall,  the  linearized  Euler  equations  support  three  downstream  propagating  modes  (As,  Xu,  and  A+) 
and  one  upstream  mode  (A_).  This  result  is  reflective  of  the  inherent  boundary  value  nature  of  the  subsonic 
Euler  equations.  The  downstream  modes  should  be  specified  at  the  domain  inlet  and  the  upstream  mode 
should  be  specified  at  the  domain  outlet. 

If  instead  the  problem  is  solved  as  in  initial  value  problem  in  space  -  by  specifying  all  modes  at  the  inlet 
and  marching  the  solution  downstream  -  the  upstream  acoustic  mode  will  cause  instability  in  the  march. 
Mathematically,  this  instability  occurs  because  A_  takes  on  complex  values  with  negative  imaginary  part  (see 
figure  1),  causing  exponential  growth  of  the  mode.  As  mentioned  in  the  introduction,  PSE  achieves  a  stable 
downstream  march  by  numerically  damping  the  upstream  mode.  This  effectively  eliminates  the  instability, 
but  also  heavily  damps  the  downstream  acoustic  modes.  Instead,  we  will  develop  a  method  which  completely 
removes  the  upstream  acoustic  modes  from  the  Euler  equations  while  retaining  all  of  the  downstream  modes. 
This  will  result  in  a  new  one-way  Euler  operator  that  accurately  models  the  downstream  Euler  modes,  but 
unlike  the  full  Euler  operator,  can  be  stably  integrated  in  the  downstream  direction. 

II. B.  Partial  decoupling  of  dynamics 

Our  method  contains  two  steps.  First,  the  upstream  and  downstream  dynamics  contained  within  the  Euler 
equations  are  partially  decoupled.  In  what  follows,  we  drop  the  subscript  n  which  indicated  the  frequency 
mode  for  notational  brevity.  Each  mode  is  treated  independently  in  the  same  manner.  Suppose  that  there 
exists  a  linear  similarity  transformation  Tdiag  such  that  the  change  of  variables 
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(ii) 


^  =  {  }  =  Tdiagq, 

upon  substitution  into  equation  (7),  results  in  a  block  diagonalized  system  of  equations 
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where  all  eigenvalues  of  D++  have  positive  group  velocity  and  all  eigenvalues  of  D _ have  negative  group 

velocity.  Then,  ij}+  contains  all  of  the  possible  downstream  modes  and  no  upstream  modes,  while  t/>_  contains 
all  of  the  upstream  modes  and  no  downstream  modes.  The  one-way  evolution  of  this  system  is  specified 
exactly  by  the  pair  of  equations 


^i/>+  =  D++^+1  (13a) 

dx 

V>_  =  0.  (13b) 

In  general,  the  transformation  required  to  block  diagonalize  the  Euler  equations  is  unknown  analytically 
and  very  expensive  to  compute  numerically.  Suppose,  instead,  that  we  choose  a  different  transformation  T 
that  does  not  diagonalize  the  system,  but  that  we  expect  to  be  a  reasonable  approximation  of  the  diagonal¬ 
izing  transformation  T^iag-  Then  inserting  the  change  of  variables 
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=  Tq , 


into  equation  (7)  will  give  an  equation  of  the  form 


with 


P++  P+- 
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(14) 
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The  off-diagonal  terms  P_| _ and  P _ j.  are  non-zero  and  nothing  is  implied  about  the  eigenstructure  of  the 

diagonal  matrices  P++  and  P _ Now,  (f>+  will  consist  mostly  of  downstream  modes,  but  will  also  contain 

some  upstream  content.  Likewise,  </>_  will  be  primarily  made  up  of  upstream  modes,  but  will  also  contain 
some  downstream  content.  In  other  words,  the  transformation  has  partially  decoupled  the  upstream  and 
downstream  dynamics. 

The  uniform  flow  Euler  eigenvectors  given  by  equation  (9)  are  a  useful  guide  for  choosing  a  transformation 
for  the  Euler  equations.  If  in  fact  we  are  working  with  a  truly  uniform  flow,  then  the  choice  T  =  F~1{U(z )}, 
with  the  function  fi(z)  replaced  by  a  convergent  rational  approximation  of  the  square  root,  will  completely 
diagonalize  the  system.  The  operator  F~l  represents  the  inverse  transverse  Fourier  transform.  This  is  the 
approach  used  in  developing  one-way  approximations  of  the  wave  equations.18, 19  For  a  non-uniform  flow,  in 
which  case  the  exact  diagonalizing  transformation  is  not  known,  the  choice  T  =  F~l {U (z)}  still  gives  a  good, 
although  incomplete,  decoupling.  The  downside  of  this  choice  is  that  the  inverse  Fourier  transform  of  the 
rational  approximation  of  the  non-local  function  fi{z)  results  in  a  P  which  is  a  full  matrix  upon  discretization 
in  y.  This  turns  out  to  be  very  undesirable  for  the  second  step  of  our  method. 

Instead,  we  use  the  even  simpler  transformation  T  =  U{z  =  0).  Under  this  transformation,  becomes 
the  well-known  one-dimensional  characteristic  variables: 
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The  characteristic  variables  are  an  entropy-like  variable  s  and  partially  decoupled  right 
like  variables  f  and  l.  The  first  three  make  up  the  approximate  upstream  variable  <fi+ 


and  left  acoustic- 
and  the  last  one 
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is  the  approximate  downstream  variable  </>_.  Their  name  comes  from  the  fact  that  this  transformation 
completely  decouples  the  upstream  and  downstream  modes  for  the  one-dimensional  Euler  equations.  In  a  two 
dimensional  setting,  the  transformation  will  perfectly  decouple  waves  traveling  directly  up  or  downstream, 
and  the  quality  will  decrease  for  waves  traveling  at  increasing  angels  from  the  axis.  This  choice  maintains 
the  sparsity  of  the  chosen  transverse  discretization. 

Since  the  upstream  and  downstream  dynamics  are  now  partially  decoupled,  a  possible  approach  would 
be  to  ignore  the  remaining  coupling  between  the  upstream  and  downstream  dynamics  and  set  <j>_  =  0.  The 
nature  of  the  error  caused  by  this  approximation  is  revealed  if  we  examine  the  relationship  between  the  true 
upstream  and  downstream  variables  and  the  approximate  ones.  Since  both  the  (hypothetical)  diagonalizing 
transformation  T^iag  and  the  approximate  diagonalization  T  are  linear  similarity  transformations,  there 
exists  a  linear  relationship  between  the  fully  decoupled  and  partially  decoupled  variables: 


</>+ 

4>- 


U+-  1  f  0V  1 

u-  \  l  +-  S' 


(18) 


We  make  two  important  observations  about  equation  (18).  First,  if  </>_  is  set  to  zero,  the  true  upstream 
and  downstream  variables  will  be  incorrect.  In  particular,  certain  upstream  solution  will  be  allowed  since 

xjj-  ^  0  and  certain  downstream  solutions  will  be  disallowed  since  the  U. \ _ </>_  contribution  to  is  ignored. 

Second,  notice  that  there  exists  a  relationship  between  the  approximate  upstream  and  downstream  variables 
4>+  and  (/>_  that  results  in  the  desired  outcome  =  0.  Precisely,  this  is  accomplished  when 


0-  =  -UZ^U-+4>+  :=  F(f>+. 


(19) 


When  this  condition  is  met,  <p  contains  only  downstream  content.  Substituting  this  result  into  equa¬ 
tion  (15)  gives: 


—  </>+  =  (P++  +  P+.F)  <t>+, 
4>-  =  F(f>  +  . 


(20a) 

(20b) 


This  result  is  an  exact  parabolization  of  the  equations  if  F  is  known  exactly.  Our  approximate  one-way 
equation  is  produced  by  finding  an  approximate  F. 


II. C.  Recursive  filtering 

To  this  end,  we  introduce  a  set  of  auxiliary  variables  :  j  =  0, . . . ,  Np}  that  satisfy  the  recursion  relation 

(■ P  -  bjlW  =(P~  h!W+\  J  =  Nf  -  1  (21) 

with  the  two-sided  termination  condition 


0+  =  </>+  (22) 

4>"F  =  0.  (23) 

Similar  recursions  were  introduced  by  Givoli  &  Neta20  and  Hagstrom  &  Warburton21  within  the  context 
of  radiation  boundary  conditions.  In  that  application,  the  goal  is  to  specify  a  boundary  condition  that  is 
transparent  to  outgoing  waves  and  does  not  generate  spurious  incoming  waves.  This  recursion  accomplish 
this  by  generating  an  operator  at  the  boundary  that  is  satisfied  by  waves  with  outgoing  group  velocities  but 
not  satisfied  by  waves  with  incoming  group  velocities.  This  is  very  similar  to  our  goal.  We  wish  to  develop 
an  operator  that  only  admits  downstream  traveling  solutions.  Our  use  of  these  recursions  differs  only  in  that 
we  use  them  to  find  the  proper  relationship  between  <p+  and  <^>_ ,  while  boundary  condition  applications  use 
them  to  relate  time  derivatives  of  these  quantities.  Additional  details  about  this  recursive  method  within 
the  context  of  boundary  conditions  can  be  found  in  the  literature. 20-24 

The  recursion  is  completely  defined  by  its  order  Np  and  the  sets  of  complex  integers  {bj .  b3  :  j  = 
0,...,Np  —  1},  which  we  call  filtering  parameters.  The  system  of  equations  defined  by  equation  (21), 
equation  (22),  and  equation  (23)  can  be  solved  for  as  a  function  of  <p+ .  If  the  filtering  parameters  are 
properly  selected,  the  recursions  converge  to  give  the  correct  relationship  between  these  variables  as  defined 
by  equation  (19).  To  achieve  convergence,  the  bj' s  and  bj' s  must  be  chosen  such  that  they  are  located  near 
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the  modes  we  wish  to  retain  and  remove,  respectively.  For  the  Euler  equations,  we  choose  these  locations 
based  on  our  knowledge  of  the  uniform  flow  eigenvalues: 

bj  £  {*As,*Aomega,*A+}  for  z  e  (0,  oo),  (24a) 

bj  £  lA _  for  z  £  (0,  oo).  (24b) 

To  account  for  baseflow  non-uniformity,  in  addition  to  varying  z,  we  also  evaluate  the  uniform  flow 
eigenvalues  using  the  range  of  k  and  Mx  that  exist  in  the  baseflow  at  the  current  x  location  in  the  march. 
This  ensures  that  the  total  set  of  filtering  parameters  will  contain  entries  that  target  all  of  the  actual  non- 
uniform  modes. 

We  will  show  in  the  next  section  that  selecting  the  filtering  parameters  based  on  equation  (24)  is  sufficient 
to  achieve  convergence.  However,  we  are  doubtful  that  this  choice  is  optimal.  Optimal  parameters  have  been 
derived  for  the  wave  equation  within  the  context  of  nonreflecting  boundary  conditions  that  guarantee  spectral 
convergence.22  These  values  are  fundamentally  linked  to  a  time-domain  formulation  and  cannot  be  easily 
translated  to  the  frequency  domain.  Our  selection  strategy  is  an  ad  hoc  frequency  domain  adaption  of  these 
parameters,  though  we  do  typically  observe  spectral  convergence.  However,  it  is  not  clear  that  optimal  time 
domain  parameters  are  optimal  in  the  frequency  domain,  and  we  have  reason  to  expect  that  they  are  not. 
Finding  optimal  frequency  domain  parameters  is  an  ongoing  aspect  of  our  work. 

II.  D.  Computational  Cost 

Implementation  of  the  recursion  equations  requires  the  solution  of  a  ANyNp  x  ANyNp  system  of  equations, 
where  Ny  is  the  number  of  discretization  points  in  y  and  Np  is  the  order  of  the  recursion.  The  number 
four  appears  because  there  are  four  dependent  variables  in  the  two-dimensional  Euler  equations.  We  use  the 
software  package  MUMPS,  which  is  a  sparse,  multi-frontal,  direct  solver. 

We  have  found  that  the  variables  in  the  system  of  equations  can  be  reordered  to  create  a  matrix  whose 
bandwidth  scales  linearly  with  Np  and  does  not  depend  on  Ny.  Therefore,  we  expect  solving  the  matrix 
system  to  require  0{NyNp)  operations.  However,  we  have  observed  that  the  actual  cost  scales  like  0(NyNp). 
In  comparison,  performing  one  PSE  step  requires  0(NyNn)  operations,  where  Nlt  is  the  number  of  iterations 
required  to  satisfy  a  nonlinear  constraint  that  is  part  of  the  PSE  formulation. 

III.  Results 

In  this  section,  we  present  results  from  two  test  problems  that  demonstrate  the  convergence  and  utility 
of  our  one-way  Euler  method.  First,  we  model  the  propagation  of  acoustic  waves  in  a  uniform  flow  created 
by  a  monopole  located  upstream  of  the  computational  domain.  Second,  we  simulate  the  generation  and 
propagation  of  acoustic  waves  in  a  convectively  unstable  turbulent  mixing  layer.  For  both  problems,  the 
one-way  Euler  results  are  compared  to  corresponding  PSE  solutions. 

For  all  simulations,  the  equations  are  discretized  in  the  transverse  direction  using  fourth-order  central 
finite  differences  and  the  computational  domain  is  truncated  using  the  super-grid  damping  layer  method25,26 
to  enforce  radiation  boundary  conditions.  For  the  one-way  Euler  equations,  we  performed  the  downstream 
march  using  fourth-order  Runge-Kutta  integration,  while  PSE  uses  implicit  Euler  integration  by  definition. 

III.  A.  Monopole  forcing 

In  this  problem,  a  monopole  disturbance  is  placed  just  upstream  of  the  computational  domain  in  a  Mach 
0.5  uniform  stream.  The  monopole  generates  an  acoustic  field  for  which  an  analytical  solution  exists.  Along 
the  inlet  of  the  computational  domain,  the  exact  solution  is  supplied  to  the  one-way  Euler  equations  and 
to  the  PSE  algorithm  as  an  initial  condition.  Then,  this  solution  is  propagated  downstream  by  the  spatial 
marching  routine  of  each  method.  Figure  2(a)  shows  the  real  part  of  the  pressure  field  created  by  the 
monopole.  The  wavefronts  are  not  concentric  because  of  the  non-zero  baseflow  moving  from  left  to  right. 
The  initial  condition  is  specified  at  the  location  of  the  dashed  line,  and  the  solution  to  the  right  of  this  is 
reproduced  by  the  marching  procedures. 

The  error  of  the  overall  full-field  solutions  produced  by  PSE  and  the  one-way  Euler  equations  are  calcu¬ 
lated  by  comparison  to  the  full-field  exact  analytical  solution.  Figure  2(b)  contains  the  error  results.  The 
convergence  of  the  one-way  Euler  solution  is  roughly  spectral  as  the  order  of  the  recursion  increases  until  the 
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error  levels  of  at  Np  =  4.  By  this  point,  the  error  in  the  one-way  Euler  solution  is  two  orders  of  magnitude 
lower  than  the  error  in  the  PSE  solution.  The  plateau  in  the  one-way  Euler  convergence  occurs  when  the 
discretization  error  becomes  larger  than  the  parabolization  error. 


Figure  2:  Monopole  test  case,  (a)  real  part  of  the  pressure  field  for  the  exact  monopole  solution.  The  dashed 
line  shows  the  location  where  the  initial  condition  is  specified  for  the  marching  procedures,  (b)  global  error 
of  one-way  Euler  equations  (symbols)  and  PSE  (line)  compared  to  the  exact  solution.  Error  is  defined  as 
||  (p  —  v exact) li> exact  H2  integrated  over  the  full  computational  domain. 


A  useful  way  to  visualize  the  convergence  of  the  parabolization  procedure  is  to  compare  the  eigenvalues 
of  the  one-way  Euler  operator  to  those  of  the  original  Euler  equations.  Precisely,  we  compare  the  eigenvalues 
of  P  (see  equation  (15))  against  those  of  the  one-way  operator  given  in  equation  (20a)  for  different  values  of 
the  order  of  the  recursion.  Figure  3  shows  these  eigenvalues  for  several  of  the  operators  used  in  the  monopole 
example  problem  above.  Note  that  additional  stable  convective  modes  exist  as  in  figure  1,  but  here  we  have 
omitted  these  eigenvalues  from  the  plot  for  readability.  Also,  the  moderate  change  in  shape  of  the  acoustic 
branches  compared  to  figure  1  is  caused  by  the  radiation  boundary  conditions,  which  are  not  enforced  for 
the  idealized  plane  waves  used  to  derive  the  analytical  modes. 

Recall  that  the  goal  of  our  method  is  to  generate  an  operator  that  accurately  contains  the  downstream 
Euler  modes  but  does  not  include  any  upstream  acoustic  modes.  That  this  goal  is  achieved  is  clearly 
demonstrated  by  figure  3.  The  Np  =  0  case  corresponds  to  performing  no  filtering  and  depending  entirely 
upon  the  partial  decoupling  of  the  dynamics  realized  by  the  transformation  to  one-dimensional  characteristic 
variables  .  Even  with  no  filtering,  the  upstream  acoustic  waves  have  been  removed.  The  very  small  angle 
downstream  propagating  waves  are  properly  captured  while  all  other  waves  are  poorly  modeled.  This  is 
consistent  with  our  previous  discussion  about  the  nature  of  the  one-dimensional  characteristic  variables. 

As  Np  is  increased,  more  and  more  of  the  downstream  acoustic  modes  are  properly  modeled.  At  Np  =  4, 
all  propagating  waves  are  well  resolved.  Increasing  Np  beyond  this  point  primarily  effects  less  important 
evanescent  modes,  although  close  examination  revels  that  all  modes  continue  to  converge  toward  the  spectrum 
of  the  original  discretized  Euler  equations.  Note  that  convergence  to  the  spectrum  of  P  implies  that  we 
have  represented  all  waves  as  well  as  possible  for  the  given  transverse  discretization.  The  properties  of  the 
underlying  discretization  are  unaltered  by  the  parabolization  method.  One  manifestation  of  this  is  that  the 
eigenvalues  of  the  parabolized  operator  converge  to  the  physical  and  unphysical  modes  that  exist  in  the 
original  discretized  Euler  equations.  This  behavior  is  evident  in  figure  3. 
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(b)  Nf  =  0 


(c)  Nf  =  2 
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Figure  3:  Convergence  of  the  spectrum  of  the  one-way  operator  to  the  downstream  modes  of  the  original 
Euler  operator.  The  upstream  waves  are  completely  removed  and  the  downstream  waves  converge  rapidly. 
The  circles  are  the  original  Euler  acoustic  modes,  and  the  pluses  are  the  one-way  Euler  acoustic  modes. 
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III.B.  Turbulent  mixing  layer 

Here,  we  model  the  generation  and  propagation  of  sound  in  a  turbulent  mixing  layer.  The  baseflow  is  a 
second-order  approximation  of  a  self-similar  solution  of  the  incompressible  boundary  layer  equations  using 
Prandtl’s  eddy  viscosity  model.2'  The  Mach  number  of  the  fast  and  slow  free-streams  are  0.8  and  0.2, 
respectively.  The  shear  layer  thickness  grows  linearly,  with  the  spread  rate  chosen  to  match  experiments. 
The  length  scale  is  chosen  such  that  the  vorticity  thickness  is  one  at  the  inlet.  We  show  results  for  the 
most  unstable  frequency  at  the  inlet.  The  computational  domain  extends  150  vorticity  thicknesses  in  the 
transverse  direction  and  100  vorticity  thicknesses  in  the  axial  direction.  We  retain  such  a  large  domain  in 
order  to  clearly  show  the  far-held  acoustics.  This  domain  is  discretized  using  3000  points  in  y.  The  one-way 
Euler  solution  uses  a  step  size  of  Ax  =  0.05  while  the  PSE  solution  uses  Ax  =  1,  which  is  just  above  the  PSE 
minimum  step  size  restriction  for  this  problem.  See  the  beginning  of  section  III  for  additional  computational 
details. 

The  initial  condition  at  the  inlet  is  a  pure  Kelvin-Helmholtz  eigenfunction  computed  at  x  =  0  under  a 
locally-parallel  assumption.  Since  the  initial  condition  does  not  include  any  acoustic  disturbances,  all  sound 
in  the  flow  is  created  by  the  growth,  decay,  and  distortion  of  the  Kelvin-Helmholtz  mode.  Because  of  the 
central  role  of  this  mode,  it  is  critical  that  the  one-way  Euler  operator  includes  it  accurately.  Figure  4  shows 
that  the  eigenvalue  of  the  Kelvin-Helmholtz  mode  in  the  one-way  Euler  operator  converges  spectrally  to  the 
eigenvalue  of  the  Kelvin-Helmholtz  mode  in  the  original  Euler  operator.  The  acoustic  modes  converge  in  a 
similar  manner.  For  the  upcoming  results,  we  use  Np  =  15. 
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Figure  4:  Convergence  of  the  one-way  Euler  Kelvin-Helmholtz  eigenvalue  to  the  original  Euler  Kelvin- 
Helmholtz  eigenvalue.  Error  is  defined  as  the  magnitude  of  the  difference  between  the  one-way  Euler  eigen¬ 
value  and  the  original  Euler  eigenvalue  normalized  by  the  magnitude  of  the  original  eigenvalue. 


Figure  5  shows  the  pressure  field  as  computed  by  the  PSE  and  the  one-way  Euler  equations.  Compar¬ 
ing  the  two  solutions  reveals  two  benefits  of  the  later  method.  First,  notice  the  difference  in  the  near-field 
wavepacket  -  the  PSE  wavepacket  decays  faster  than  the  one-way  Euler  wavepacket.  These  near-field  dy¬ 
namics  are  dominated  by  the  growth  and  decay  of  the  Kelvin-Helmholtz  mode.  We  have  shown  in  figure  4 
that  the  Kelvin-Helmholtz  mode  in  the  one-way  operator  is  converged  to  the  correct  value.  At  Np  =  15,  the 
error  is  less  than  10~9,  so  the  accumulated  parabolization  error  over  the  length  of  the  domain  is  less  than 
exp(100  x  10~9)  «  10“ '.  The  error  due  to  the  discretization  is  also  small  because  of  the  small  step  size  and 
fourth-order  integration. 

The  overly  rapid  decay  of  the  near-field  wavepacket  in  the  PSE  solution  can  be  explained  by  the  implicit 
Euler  integration  of  the  PSE  equations.  In  this  flow,  the  Kelvin-Helmholtz  mode  is  most  unstable  very  near 
the  inlet,  so  the  growth  rate  is  decreasing  (or  the  decay  rate  is  increasing)  over  the  majority  of  the  domain. 
Therefore,  over  each  step  in  the  PSE  march,  the  growth  rate  is  lower  (or  the  decay  rate  is  higher)  at  the 
downstream  side  of  the  interval  than  at  the  upstream  side.  Because  of  the  implicit  Euler  integration,  the 
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(a)  (b) 

Figure  5:  Pressure  field  from  the  turbulent  mixing  layer  test  case  for  (a)  PSE  and  (b)  the  one-way  Euler 
equation. 


downstream  value  is  applied  over  the  entire  interval,  so  the  evolution  of  the  mode  is  always  biased  toward 
a  more  damped  solution.  As  the  step  size  decreases,  the  integration  converges  toward  the  correct  result, 
but  satisfactory  convergence  cannot  be  achieved  prior  to  the  violation  of  the  minimum  step  size  restriction 
required  for  stability  of  the  PSE  march. 

Second,  the  one-way  Euler  solution  contains  a  well-resolved  acoustic  field  while  the  PSE  solution  contains 
almost  no  acoustic  field  at  all.  The  downstream  acoustic  modes  in  the  PSE  equations  were  heavily  damped 
in  order  to  adequately  damp  the  upstream  acoustic  modes  to  allow  a  stable  march.  The  downstream 
acoustic  modes  in  the  one-way  Euler  operator,  on  the  other  hand,  were  properly  retained,  allowing  the 
strong  downstream  radiation  produced  by  the  peaking  and  subsequent  decay  of  the  Kelvin-Helmholtz  mode 
to  be  captured. 

In  figure  6(b),  we  show  a  slice  of  the  pressure  field  along  x  =  25.  In  the  near-field,  the  PSE  and 
one-way  Euler  results  are  similar,  with  the  minor  discrepancy  caused  by  the  growth  rate  error  discussed 
above.  The  exponential  decay  of  the  near-field  pressure  is  consistent  with  the  asymptotic  behavior  of  the 
Kelvin-Helmholtz  mode.  The  one-way  Euler  pressure  field  levels  off  due  to  the  presence  of  propagating 
(non-decaying)  acoustic  waves  that  radiate  to  the  far-held.  In  contrast,  the  damped  remnants  of  acoustic 
radiation  in  the  PSE  solution  is  at  least  four  order  of  magnitude  too  small.  These  results  are  consistent  with 
comparisons  made  between  DNS  and  PSE  by  Cheung  &  Lele28  for  laminar  subsonic  mixing  layers. 

Lastly,  we  address  the  irregular  waveforms  visible  in  the  one-way  Euler  solution  at  high  very  high  angles 
near  x  =  0.  This  error  is  due  to  slow  convergence  of  a  hand-full  of  high  angle  acoustic  modes.  These  modes 
do  converge  as  the  order  of  the  filtering  is  increased,  but  more  slowly  than  the  other  modes  in  the  Euler 
spectrum.  Therefore,  this  error  can  be  removed,  but  the  filtering  order  required  to  do  so  is  much  higher  than 
necessary  for  all  other  modes.  Eliminating  this  non-uniform  convergence  is  one  of  the  goals  of  our  ongoing 
efforts  to  derive  optimal  filtering  parameters. 
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(a)  (b) 

Figure  6:  One-way  Euler  and  PSE  pressure  amplitude  for  the  turbulent  mixing  layer  (a)  along  the  axial 
slice  y  =  0  (b)  along  the  transverse  slice  x  =  25. 


IV.  Summary 

We  have  developed  a  one-way  Euler  equation  that  can  be  efficiently  and  stably  solved  as  a  spatial  initial 
value  problem.  The  one-way  operator  contains  no  upstream  propagating  modes,  and  can  represent  all 
downstream  modes  with  arbitrary  accuracy. 

For  problems  where  the  mid-  and  far-field  acoustic  radiation  is  of  interest,  this  technique  will  provide 
significant  improvements  over  PSE.  For  example,  we  plan  to  use  this  technique  as  a  reduced  order  model  to 
obtain  efficient  jet  noise  predictions. 
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